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Abstract: We evaluate the relic density of neutralinos in the minimal supergravity 
(mSUGRA) model. All 2 — * 2 neutralino annihilation diagrams, as well as all processes 
involving sleptons, charginos, neutralinos and third generation squarks are included. Rel- 
Qlh! ativistic thermal averaging of the velocity times cross sections is performed. We find that 

co-annihilation effects are only important on the edges of the model parameter space, where 
some amount of fine-tuning is necessary to obtain a reasonable relic density. Alternatively, 



at high tan /?, annihilation through very broad Higgs resonances gives rise to an acceptable 
neutralino relic density over broad regions of parameter space where little or no fine-tuning 
is needed. Finally, we compare our results against the reach of various e + e~ and hadron 
colliders for supersymmetric matter. 



1. Introduction 



A wide variety of astrophysical measurements are being used to pin down some of the basic 
cosmological parameters of the universe. High resolution maps of the cosmic microwave 
background (CMB) radiation []IJ] imply that the energy density of the universe Q = p/p c — 1, 
consistent with inflationary cosmology. Here, p c = 3H 2 /8ttGn is the critical closure density 
of the universe, where Gn is Newton's constant and H = lOO/i km/sec/Mpc is the scaled 
Hubble constant. The value of h itself is determined to be ~ 0.7 ±0.1 by improved 
measurements of distant galaxies^. Meanwhile, data from distant supernovae||] imply 
a nonzero dark energy content of the universe ~ 0.7, a result which is confirmed by 
fits to the CMB power spectrum^]]. Analyses of Big Bang nucleosynthesis^] imply the 
baryonic density Q^h 2 ~ 0.020 ± 0.002, although the CMB fits suggest a somewhat higher 
value of ~ 0.03. Hot dark matter, for instance from massive neutrinos, should give only 
a small contribution to the total matter density of the universe. In contrast, a variety of 
data ranging from galactic rotation curves to large scale structure and the CMB imply a 
significant density of cold dark matter (CDM)|| Q c h 2 ~ 0.2 ± 0.1. 

In many i?-parity conserving supersymmetric models of particle physics, the lightest 
neutralino [Z\) is also the lightest SUSY particle (LSP); as such, it is massive, neutral and 
stable. For this case, relic neutralinos left over from the Big Bang provide an excellent 
candidate for the CDM content of the universe^. In the early universe, such neutralinos 
would exist in thermal equilibrium with the cosmic soup. As the universe expanded and 
cooled, the thermal energy would no longer be sufficient to produce neutralinos at an ap- 
preciable rate, although they could still annihilate away. Their number density is governed 
by the Boltzmann equation formulated for a Friedmann-Robertson- Walker universe. 

In this paper, our goal is to present results of calculations of the neutralino relic 
density within the context of the paradigm minimal supergravity model (mSUGRA, or 
CMSSM) [|| . In mSUGRA, it is assumed that SUSY breaking occurs in a hidden sector of 
the model, with SUSY breaking effects communicated from hidden to observable sectors 
via gravitational interactions. The model parameter space is given by 

m 0; ^1/2) A)) tan/3 and sign(p). (1-1) 

Here, mo is the universal scalar mass, m 1/ / 2 is the universal gaugino mass and Aq is the 
universal trilinear mass all evaluated at Mqut, while tan/3 is the ratio of Higgs field vevs 
(v u /vd), and p, is a supersymmetric Higgs mass term. The soft SUSY breaking parameters, 
along with gauge and Yukawa couplings, evolve from Mqut to M wea k according to their 
renormalization group (RG) equations. At M wea k, the RG improved 1-loop effective po- 
tential is minimized, and electroweak gauge symmetry is broken radiatively. In this report, 
we implement the mSUGRA solution encoded in ISAJET v7.59||. 

There is a long history of increasingly sophisticated solutions for the relic density of 
neutralinos in supersymmetric models go|, |l|, ||, ||, [l|, 0, H, HI, H, H H 



24], |25|, |26|, |27], |28|, |29], |3^, [31], |32[| . The key ingredient to solving the Boltzmann equation 
is to evaluate the thermally averaged neutralino annihilation cross section times velocity 
factor. Traditionally, the solution is made by expanding the annihilation cross section as a 
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power series in neutralino velocity, so that angular and energy integrals can be evaluated 
analytically. The remaining integral over temperature can then be performed numerically. 
The power series solution is excellent in many regions of model parameter space because 
the relic neutralino velocity is expected to be highly non-relativistic. 

However, it was emphasized by Griest and Seckel that annihilations may occur through 
s-channel resonances at high enough energies [14] that a relativistic treatment of thermal 
averaging might be necessary. Drees and Nojiri found that at large values of the parameter 
tan P, neutralino annihilation can be dominated by s-channel scattering through very broad 
A and H Higgs resonances [17]. The proper formalism for relativistic thermal averaging 
was developed by Gondolo and Gelmini (GG) |i~5|| , and was implemented in the code of 



Baer and BrhlikpO], 22]. Working within the framework of the mSUGRA model, it was 



found[2C, 22, 23, |2J|, |3Cfl that at large tan/?, indeed large new regions of model parameter 
space gave rise to reasonable values for the CDM relic density. At large tan j3, the A and H 
resonances are broad enough (typically 10-50 GeV) that even if the quantity 2m^ i is several 
partial widths away from exact resonance, there can still be a significant rate for neutralino 
annihilation. Thus, in the mSUGRA model at low mo and tan /3, neutralino annihilation is 
dominated by i-channel slepton exchange, and reasonable values of the relic density occur 
only for relatively low values of tuq and mi/2- At high tan/?, a much larger parameter 
space is allowed, owing to off-resonance neutralino annihilation through the broad Higgs 
resonances. 

In addition, there exist regions of mSUGRA model parameter space where co-an- 
nihilation processes are important, and even dominant. It was stressed by Griest and 
Seckel[14] that in regions with a higgsino-like LSP, the Z±, W\ and Z<i masses become 
nearly degenerate, so that all three species can exist in thermal equilibrium, and annihilate 
against one another. The relativistic thermal averaging formalism of GG was extended to 
include co-annihilation processes by Edsjo and Gondolo |2"lf|, and was implemented in the 
DarkSUSY code] 24] for co-annihilation of charginos and heavier neutralinos. 

The importance of neutralino-slepton co-annihilation was stressed by Ellis et al. and 
others(2^, 27, |28|, p0| . In regions of mSUGRA parameter space where Z\ and f\ (or 
other sleptons) were nearly degenerate (at low mo), co-annihilations could give rise to 
reasonable values of the relic density even at very large values of m^, at both low and 
high tan/3. In addition, for large values of the parameter Aq or for non-universal scalar 
masses, top or bottom squark masses could become nearly degenerate with the Z\, so that 
squark co-annihilation processes can become important as well|3l], 32]. 

In this paper, our goal is to calculate the relic density of neutralinos in the mSUGRA 
model including co-annihilation processes in addition to relativistic thermal averaging of 
the annihilation cross section times velocity. Since there are very many Feynman diagrams 
to evaluate for neutralino annihilations and co-annihilations, we use CompHEP v. 33.24 
which provides for fast and efficient automatic evaluation of tree level processes in the SM 
or MSSM. For initial states including Z\, Z2, W\, e\, fii, f%, i\ and 61, we count 1722 
subprocesses, including 7618 Feynman diagrams. For those processes we have calculated 
the squared matrix element and have written it down in the form of CompHEP FORTRAN 
output. 
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The weak scale parameters from super symmetric models are generated using ISAJET 
v7.59, and interfaced with the squared matrix elements from CompHEP. Details of our 
computational algorithm are given in Sec. 2. In Sec. 3 we present a variety of results for 
the relic density in mSUGRA model parameter space. Much of parameter space is ruled 
out at low tan j3 since the relic density is too high, and would yield too small an age of the 
universe. At high tan /?, large regions of parameter space are available with a reasonable 
relic density in the range 0.1 < % h 2 < 0.3. In Sec. 4, we compare our results with 
some previous results on the reach of colliders, and draw some implications. In Sec. 5, we 
conclude. 

As we were completing this work, the group of Belanger, Boudjema, Pukhov and 



Semenov reported on a calculation similar to ours in scope and method[34[. In addition 



a paper by Nihei, Roszkowski and de Austri appeared, containing analytic calculations of 



all ZiZi annihilation cross sections [35]. 
2. Calculational Details 

The evolution of the number density of super symmetric relics in the universe is described 
by the Boltzmann equation as formulated for a Friedmann-Robertson- Walker universe. For 
calculations including many particle species, such as the case where co-annihilations are 
important, there is a Boltzmann equation for each particle species. Following Griest and 
Seckel[|l4|, the equations can be combined to obtain a single equation 

^ = _ 3Hn _ {aeffv) ( n a _ n y ( 2.i) 

where 

N 

n = ^2m (2.2) 

i=l 

and the sum extends over the N particle species contributing to the relic density, with m 
being the number density of the ith species. Furthermore, n eq j is the number density of 
the ith species in thermal equilibrium, given by 

9im 2 T /rm\ , ■ 

ne q ,i = ^ r K 2 l K —) , (2.3) 

where Kj is a modified Bessel function of the second kind of order j. 

The quantity (a e ffv) is the thermally averaged cross section times velocity. A succinct 
expression for this quantity using relativistic thermal averaging was computed by Gondolo 



and Gelmini for the case of a single particle species [15], and was extended by Edsjo and 



Gondolo for the case including co- annihilations p, 1| . We adopt this latter form, given by 

_ r*(»x^v.»»»y»w» (2 . 4) 



^(e£,a' 2 (|)% 



where x = T/m^^ is the temperature in units of mass of the relic neutralino, Ojj is the cross 
section for the annihilation reaction ij — > X (X is any allowed final state consisting of 2 
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SM and/or Higgs particles), A(a 2 , b 2 , b 2 ) = a 4 + bf + tf- 2{a 2 b 2 + a 2 b 2 + b 2 b 2 ), a = y/s/m^ 
and bi = mj/m^. This expression is our master formula for the relativistically thermal 
averaged annihilation cross section times velocity. 

To solve the Boltzmann equation, we introduce a freeze-out temperature Tp, so that 
the relic density of neutralinos is given by 1 

fi h * = PP°1 ( 2 .5) 

Zl 8.1 x 10- 47 GeV 4 v ' 

where 

M Pi y T 7 y J Q (a e ffv)dx 

The freeze-out temperature xp = Tp/m^ is determined as usual by an iterative solution 
of the freeze-out relation 



log 



m z 1 9eff I 45 , ,1/2 



(2.7) 



2^3 2 \l2g if G N V" t »"'^ 1,,J ' F 
Here, g e ft denotes the effective number of degrees of freedom of the co-annihilating par- 



ticles, as defined by Griest and Seckel[14|. The quantity g* is the SM effective degrees of 
freedom parameter with ^fgl ~ 9 over our region of interest. 

The challenge then is to evaluate all possible channels for neutralino annihilation to SM 
and/or Higgs particles, and also all co-annihilation reactions. The 7618 Feynman diagrams 
are evaluated using CompHEP, leading to about 50 MB of FORTRAN code. To achieve 
our final result with relativistic thermal averaging, a three-dimensional integral must be 
performed over i.) the final state subprocess scattering angle 9, ii.) the subprocess energy 
parameter a = i/s/m^ , and Hi.) the temperature T from freeze-out Tp to the present 
day temperature of the universe, which can effectively be taken to be 0. We perform the 
three-dimensional integral using the BASES algorithm j3(|, which implements sequentially 
improved sampling in multi-dimensional Monte Carlo integration, generally with good 
convergence properties. We note that the three-dimensional integration appearing in the 
case of our relativistic calculations involving several species in thermal equilibrium is about 
2 orders of magnitude more CPU-time consuming than the series expansion approach, 
which requires just one numerical integration. 



3. Results 

Our first results in Fig. |l| show regions of Q^^ 2 in the mo vs. m 1 / 2 plane in the minimal 
supergravity model for = 0, tan (3 = 10 and for [i < and fi > 0. The upper plots show 
the contribution if only Z\Zi annihilation reactions occur, while the lower frames include 
as well all co-annihilation processes. The red shaded regions are excluded by theoretical 
constraints (lack of REWSB on the right, a charged LSP in the upper left). The unshaded 

1 The procedure we follow gives numerical results valid to about 10% versus a direct numerical solution 
of the Boltzmann equation Jl5|. 
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regions have h? > 1, and should be excluded, as they would lead to a universe of age 
less than 10 billion years, in conflict with the oldest stars found in globular clusters. The 
light blue shaded regions have ^^^i 2 < 0.02, which wouldn't be enough CDM even to 
explain galactic rotation curves. The green region yields values of 0.1 < ^^Ji 2 < 0.3 



i.e. 



in the most cosmologically favored region. The yellow (0.02 < ft^/i 2 < 0.1) and dark 
blue (0.3 < Qz_h 2 < 1) correspond to regions with intermediate values of low and high 

relic density, respectively. Points with m^ 2 ~ 150 GeV give rise to chargino masses below 
bounds from LEP2; the LEP2 excluded regions due to chargino, slepton and Higgs searches 
are not shown on these plots, but will be shown in Sec. 4. 

The structure of these plots can be understood by examining the thermally averaged 
cross section times velocity, integrated from zero temperature to Tp . In Fig. [2] we show this 
quantity for a variety of contributing subprocesses plotted versus mo for fixed m 1 / 2 = 300 
GeV, fi > 0, and all other parameters as in Fig. |l|. At low values of mo, the neutralino 
annihilation cross section is dominated by i-channel scattering into leptons pairs, as shown 
by the black solid curve. However, at the very lowest values of mo, the annihilation rate 
is sharply increased by neutralino-stau and stau-stau co-annihilations, leading to very low 
relic densities where ~ mf 1 [26|. As mo increases, the slepton masses also increase, 
which suppresses the annihilation cross section, and the relic density rises to values ^^^ 2 > 
1. When mo increases further, to beyond the ~ 1 TeV level, and approaches the excluded 
region, the magnitude of the [i parameter falls, and the higgsino component of Z\ increases. 
This is the so called "focus point" region, explored in Ref. [25]. In this region, the 
annihilation rate is dominated by scattering into WW, ZZ and Zh channels. At even 
higher mo values, ~ m w ~ m z 2 ' an< ^ these co-annihilation channels increase even 
more the annihilation rate. Finally, at the large mo bound on parameter space, |/i| — > 0, and 
appropriate REWSB no longer occurs. Most of the structure of Fig. [l] can be understood 
in these terms, with the exception being the horizontal band of very low relic density at 
m 1/2 ~ 125 GeV. 

In this region, which is nearly excluded by LEP2 bounds on the chargino mass, there 
is enhanced neutralino annihilation through the Z and h resonances. In fact, a higher 
degree of resolution on our plots would resolve these horizontal bands into two bands, 
corresponding to each of the separate resonances, as shown in Ref. p0|| . Finally, the glitch 
in contours around mo ~ 1500 GeV and m^j ~ 425 GeV occurs because — m% = 175 
GeV, so that a(ZiZi — > tt) becomes large. 

The mo vs. m-i/2 planes for tan/3 = 30 are shown in Fig. ||. The structure of these 
plots are qualitatively the same as in Fig. |l[ Quantitatively, they differ mainly in that the 
cosmologically favored regions are expanding as tan (3 grows. One reason is that the light 
stau becomes even lighter as tan j3 increases, and this increases the neutralino annihilation 
rate Z\Z\ —* rf through t-channel stau exchange. In addition, the bottom and tau Yukawa 
couplings increase with tan /3, which increases the annihilation cross sections into rs and 6s. 
Finally, the H and A Higgs boson masses are decreasing with tan (3, and annihilation rates 
which proceed through these resonances increase. Co-annihilations again gives enhanced 
annihilation cross sections on the left and farthest right hand sides of the allowed parameter 
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Figure 1: Regions of neutralino relic density in the ttiq vs. m\/2 plane for A$ = and tan/3 = 10. 
The upper two frames show the contribution for only Z\Z\ annihilation, while the lower frames 
include as well all co-annihilation processes. 



space. 

In Fig. we show the mo vs. mi/2 plane for tan/? = 45. In this case, the structure 
of the plane is changing qualitatively, especially for [i < 0. First, there is a new region 
of disallowed parameter space for fi < in the lower left due to m 2 A < 0, which signals a 
breakdown of the REWSB mechanism. Second, a corridor of very low relic density passes 
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Figure 2: Thermally averaged cross section times velocity integrated from T — to Tp, for 
various component subprocess cross sections. The blue curve denotes the total of all annihilation 
and co-annihilation reactions. We show results versus m for m X j 2 = 300 GcV, fi > 0, A = and 
tan/5 = 10. 



diagonally through the plot. The center of this region is where 2m^ — tua and ran- At 
the A and H resonance, there is very efficient neutralino annihilation into bb final states. 
This is illustrated in Fig. ^, where we show the integrated annihilation cross section times 
velocity versus mo for = 600 GeV and fj, < 0. At the very lowest values of mo, 

there is again the sharp peak due to neutralino-stau and stau-stau co-annihilations. For 
larger values of mo, however, the annihilation rate is dominantly into bb final states over 
almost the entire mo range. This is due to the large annihilation rates through the s- 
channel A and H diagrams, even when the reactions occur off resonance. In this case, the 
widths of the A and H are so large (both ~ 10 — 40 GeV across the range in mo shown) 
that efficient s-channel annihilation can occur throughout the bulk of parameter space, 
even when the resonance condition is not exactly fulfilled. The resonance annihilation is 
explicitly displayed in this plot as the annihilation bump at mo just below 1300 GeV. 
Another annihilation possibility is that Z\Z\ — > bb via t and u channel graphs. In fact, 
these annihilation graphs are enhanced due to the large b Yukawa coupling and decreasing 
value of m^ , but we have checked that the s-channel annihilation is still far the dominant 
channel. Annihilation into rr is the next most likely channel, but is always below the level 
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Figure 3: Regions of neutralino relic density in the mo vs. m\/% plane for Aq = and tan/3 = 30. 
The upper two frames show the contribution for only Z^Zi annihilation, while the lower frames 
include as well all co-annihilation processes. 

of annihilation into bb for the parameters shown in Fig. || At even higher values of uiq 
where the higgsino component of Z\ becomes non-negligible, the annihilations into WW 
and ZZ again dominate; finally, at the highest values of mo, the W\ and Z2 co-annihilation 
channels become important. 

In Fig. ||, we show again the subprocess annihilation rates versus ran for tan f3 = 45, 



-8- 



•0.3 <Qh?< 1. 0.02 <Qh*<0.1 

0.1 <Qh 2 < 0.3 a h 2 < 0.02 




oj 1000 2000 3000 1000 2000 3000 




1000 2000 3000 1000 2000 3000 

m <GeV) 



Figure 4: Regions of neutralino relic density in the mo vs. m\i2 plane for A$ = and tan/3 = 45. 
The upper two frames show the contribution for only Z±Zi annihilation, while the lower frames 
include as well all co-annihilation processes. 

but this time for /i > and for m\j2 = 300 GeV. Although no explicit resonance is evident 
for jj, > 0, the dominant annihilations are once again into bb final states over most of the 
parameter space, due to the very wide Higgs resonances. At the highest values of mo, 
where [i is becoming small, the annihilation rate into the dominant WW and ZZ final 
states becomes suppressed. The suppression is due to the diminishing mass of the Z\ as 
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Figure 5: Thermally averaged cross section times velocity evaluated at Tp for various component 
subprocess cross sections. The blue curve denotes the total of all annihilation and co-annihilation 
reactions. We show results versus mo for m^ 2 = 600 GeV, \i < 0, Aq = and tan/3 = 45. 

fj, — ► 0. As falls below Mz and then Mjy, there is thermal suppression of annihilation 
into the ZZ and WW final states. 

To summarize the regions of mSUGRA model parameter space with reasonable values 
of neutralino relic density, we can label four important regions: i.) annihilation through 
t-channel slepton- especially stau- exchange, as occurs for low values of tuq and mi/ 2 , u.) 
the stau co-annihilation region for low values of ijiq on the edge of the excluded region, 
in.) the large mo region with non-negligible higgsino-component annihilation, and also W\ 
(and possibly Z%) co-annihilation occurs near the edge of the limit of parameter space, 
and iv.) annihilation into bb and rf final states through s-channel A and H resonances at 
high tan (3. Other regions can include top or bottom squark co-annihilation for large values 
of Aq, again on the edge of parameter space where t\ or b\ become light, or annihilation 
through Z or h resonances. These latter regions are essentially excluded now by constraints 
on sparticle masses from LEP2. 

It is useful to view the relic density ^j?^ 2 directly as a function of model parameters. 
We show in Fig. ^ the value of ^^ x ^ 2 versus the parameter uiq for fixed my 2 = 600 GeV, 
Aq = 0, /i < and for tan/3 = 10, 30 and 45. The dashed curves show the result with no 
co-annihilations, while the solid curves yield the complete calculation. The shaded band 
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Figure 6: Thermally averaged cross section times velocity evaluated at Tp for various component 
subprocess cross sections. The blue curve denotes the total of all annihilation and co-annihilation 
reactions. We show results versus mo for m^ 2 = 300 GeV, \i > 0, A = and tan/3 = 45. 

denotes the cosmologically favored region with 0.1 < f^/i 2 < 0.3. For this value of m\/i, 
the lower tan (5 curves yield a favored relic density only in the very low and very high mo 
regions, and here the curves have a very sharp slope. The large slope is indicative of large 
fine-tuning, in that a small change of model parameters, in this case m-o, yields a large 
change in fl^^h 2 . In contrast, the tan/? = 45 curve shows a large region with good relic 
density and nearly zero slope, hence very little fine-tuning. 

In Fig. [8|, we show the corresponding values of the fine-tuning, basically the logarithmic 



derivative, as advocated by Ellis and Olive [37]: 

A(m ) 



mo Z\ 



h? drriQ 



(3.1) 



As indicated earlier, the low fine-tuning regions mostly coincide with that of neutralino an- 
nihilation via i-channel slepton exchange (region i.)), or off-resonance annihilation through 
A and H (region iv.)). The co-annihilation region ii.) and focus point region Hi.) tend 
to have higher fine-tunings due to the steep rise of the cross sections. Regions with si- 
multaneous low fine- tuning and preferred h 2 values are the best candidates for viable 
mSUGRA parameters. 



- 11 - 




250 500 750 1000 1250 1500 1750 2000 2250 2500 

m (GeV) 

Figure 7: Neutralino relic density fij^/i 2 versus mo for mi/a = 600, Aq — 0, \i < and tan (3 = 10, 
30 and 45. The shaded region denotes the cosmologically favored values of Q.^.Ji 2 . 

In Fig. H|, we show Q^h 2 versus itiq for = 300 GeV, ^4o = 5 M > an d the 

same three tan (3 parameters. The curves reflect the broad regions of parameter space with 
reasonable relic density values at high tan (3. The corresponding plot of the fine-tuning 
parameter is shown in Fig. [[(]. Again, there is large fine-tuning at the edges of parameter 
space, but low fine-tuning in the intermediate regions. In conclusion, the relic density and 
the fine-tuning parameter together tend to prefer mSUGRA model parameters in regions 
i.) or iv.). These two regions lead to distinct collider signatures for future searches for 
supersymmetric matter. 

4. Comparison with collider reaches 

It is worthwhile to compare our results on the neutralino relic density with various collider 
reaches 2 . To do so, we first show in Fig. |ll| the tuq vs. plane for tan/3 = 10, Aq = 
and \x > 0, but this time including information from different collider projections. First, 
the region excluded by LEP2 sparticle searches is shown by the pink shading, and reflects 
mSUGRA model points where < 100 GeV, m gl < 100 GeV, or m fl < 76 GeV|§. 

2 Non-accelerator direct and indirect dark matter search reach results are summarized in the last two 
papers of Ref. ||] 
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Figure 8: The fine tuning parameter as the function of mo for tan/3 = 10,30,45 and for the 
parameter slice mi/ 2 = 600 GeV, \i < 0. 



These LEP2 bounds sharply constrain the regions where neutralino annihilation could 
occur via the Z and light Higgs h resonances. In addition, we plot contours of light Higgs 
boson mass nih = 110, 115 and 120 GeV. Since the light Higgs scalar h is usually SM-like 
in the mSUGRA model, the region below rrih = 115 GeV is largely excluded[39] by the 
direct LEP2 Higgs search. We note that the Higgs mass varies slowly in parameter space, 
so a small change in Higgs mass can lead to large changes in model parameters. Thus, 
these bounds may have some fuzziness to them, reflecting uncertainties on the theoretical 
calculations and experimental search results. The reach of the Fermilab Tevatron for SUSY 
particles has been estimated recently in Ref. [40] for an integrated luminosity of 25 fb _1 . 
Almost all the reach comes from the search for clean trilepton events. The 3<r reach is 
denoted by the two lower black contours. The results show that Tevatron experiments will 
be able to probe a significant part of the favored relic density region where annihilation 
occurs through i-channel slepton exchange. Also, some of the "focus point" regional]] with 
large mo and small m\ii is accessible. 

The reach of the CERN LHC is also shown for 10 fb _1 of integrated luminosity [p2|j . 
The LHC reach extends well beyond the i-channel slepton region, but cannot exclude all 
the low and high too regions corresponding to slepton co- annihilations or to higgsino-like 
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Figure 9: Neutralino relic density S^/i 2 versus mo for Aq 
tan = 10, 30 and 45. 



0, m 1/2 = 300 GeV, fi > and 



neutralinos, where annihilation cross sections are enhanced. We remark, however, that 
these regions on the edge of parameter space, although perhaps not directly accessible to 
LHC searches, are also disfavored by fine-tuning requirements. 

We also show the reach of a linear e + e~ collider for SUSY particles for yfs = 500 GeV 
(NLC500) and ^ = 1000 GeV (NLC1000), assuming 30 fbr 1 of integrated luminosity @. 
The left-most NLC region is explorable via slepton pair searches, while the lower and right 
NLC regions are explorable via chargino pair searches. A small intermediate region is 
accessible via e + e~ — ► Z1Z2 searches. 

A similar comparison of neutralino relic density versus collider searches is shown in 



Fig. 12, although in this case, results for the NLC reach are unavailable. A large part of 



the green region is actually excluded by the LEP2 Higgs search results. Furthermore, the 
reach of the Fermilab Tevatron barely extends to the cosmologically favored region. The 
CERN LHC covers most of the green region, with the exception of the stau co-annihilation 
band, and the higgsino-like LSP band. 

5. Conclusions 



In conclusion, we have performed a calculation of the neutralino relic density in the min- 
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Figure 10: The fine tuning parameter as the function of mo for tan/3 = 10,30,45 and for the 
parameter slice mi/ 2 = 300 GeV, [i > 0. 



imal supergravity model including all 2 — > 2 neutralino annihilation and co-annihilation 
processes, where the initial state includes Z\, Zi-, W\, e%, fli, fi, t\ and b\. The calcu- 
lation was performed using the CompHEP program for automatic evaluation of Feynman 
diagrams, coupled with ISAJET for sparticle mass evaluation in the mSUGRA model, and 
for standard and supersymmetric couplings and decay widths. We implemented relativis- 
tic thermal averaging, which is especially important for evaluating the relic density when 
resonances in the annihilation cross section are present, and neutralino thermal velocities 
can be relativistic. The three-dimensional integration was performed by Monte Carlo eval- 
uation with importance sampling, which yields in general good convergence even in the 
presence of narrow resonances. We note once again that a calculation of similar scope and 
procedure was recently reported in Ref. [34|. 

It may be useful to compare our results against other recent evaluations of the neu- 
tralino relic density including all co-annihilation effects. We compared against several 
recently published results |2(| |23|, ^9], 30 1. Our results agree qualitatively with these other 
published results. Quantitatively, these various papers largely disagree amongst themselves 
and with our work at very large values of tan (3. In this region, even one loop calculations 
of sparticle and Higgs boson masses can be very unstable, and especially the value of is 
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Figure 11: Regions of relic density in the mo vs. m l / 2 plane, including theoretical and experimental 
constraints, contours of light Higgs mass mh (red), and reach projections for the Fermilab Tevatron, 
CERN LHC and Next Linear Collider. We adopt tan/3 = 10, A = and fi > 0. 



very sensitive to the exact procedure involved in calculating sparticle masses. This results 
in differences in the precise location of the corridor of annihilation through the A and H 
resonances. Our value of at large tan (5 seems generally larger than the values obtained 

Clearly, 



in Refs. [26, 23|, and somewhat smaller than those obtained in Refs. [29] and p0[ 
more theoretical work needs to be done to stabilize the SUSY and Higgs particle mass 
predictions at large tan (3. Finally, the width of the bands of cosmologically favored relic 
density around the A, H annihilation corridors appears much wider in our results than in 
the results of Refs. [^, 26]. This might be an effect of our improved treatment including 
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Figure 12: Regions of relic density in the mo vs. m l / 2 plane, including theoretical and experimental 
constraints, contours of light Higgs mass (red), and reach projections for the Fermilab Tevatron 
and CERN LHC. We adopt tan/3 = 45, A = and fj, > 0. 



relativistic thermal averaging. It would be useful to have a comparison against similar 
results from the group Belanger et al.[34] when these become available. 

We presented all our results within the framework of the mSUGRA model. We found 
four regions of parameter space that led to relic densities in accord with results from cos- 
mological measurements, i.e. 0.1 < ^^_h 2 < 0.3. These include i.) the region dominated 
by i-channel slepton exchange, ii.) the region dominated by stau co-annihilation, Hi.) the 
large mo region dominated by a more higgsino-like neutralino and iv.) the broad regions 
at high tan/3 dominated by off-shell annihilation through the A and H Higgs boson reso- 
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nances. Regions ii.) and Hi.) generally have large fine-tuning associated with them, and 
although it is logically possible that nature has chosen such parameters, any slight devia- 
tion of model parameters would lead to either too low or too high a relic density. Region 
i.) generally has the property that some of the sleptons have masses less than about 300- 
400 GeV. This region can give rise to a rich set of collider signatures, since many of the 
sparticles are relatively light. 

Region iv.) gives broad regions of model parameter space with reasonable values of 
relic density as well as low values of the fine-tuning parameter. It can also allow quite heavy 
values of SUSY particle masses, which would be useful to suppress many flavor-violating 
(such as b — > S7)[44| and CP violating loop processes, and the muon g — 2 value |^{|. In 
many respects region iv.) is a favored region of parameter space. The neutralino relic 
density may well point the way to the sort of SUSY signatures we should expect at high 
energy collider experiments. 
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